// Sample.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <stdexcept>
#include <iostream>
#include <map>
#include <ctime>
#include <cstdlib>
#include "string.h"
#include "blinalg.hpp"
#include "superelement.hpp"
#include "mesh.hpp"
#include "gauss.hpp"
#include "Intergration.hpp"
#include "functions.hpp"
#include "Sample.h"
#include "Solver.h"
#include <cassert>
using namespace std;
//#define Max 2
#define RefineOrder 10
#define GaussDegree 20
double A, B;
int N;
void readInput()
{
	std::ifstream txt_file;
	// open files
	txt_file.open("fem.txt", std::ios::in);
	assert(txt_file.good ());
	// read mesh
	txt_file >> N >> A >> B;
	txt_file.close();
}
int _tmain(int argc, _TCHAR* argv[])
{
	 clock_t t1, t2, t3, t4,t5,t6;
	double tol = pow(10.0, -5.0);
	readInput();
	std::cout << "Solving Elliptic Problem: -(pu')'+qu=f by FEM" << std::endl;	
	vector<FElement> m_elems;
	//Khoi tao cac pthh
	khoiTaoLuoi(m_elems,"fem");
	  t1 = clock();
	lamMinLuoi(m_elems, RefineOrder);
	  t2 = clock();
	    cout << "MakeGrid consumes " << (t2-t1)/CLOCKS_PER_SEC<< endl;
	N = m_elems.size();
	EVector<double> F(N-1),R(N-1);
	std::cout<<"Number of elems: "<<N<<endl<<endl;
	FullMatrix<double> StiffnessMatrix(N-1, N-1); //MassMatrix(N-1,N-1), 
	//Tinh ma tran B
	StiffnessMatrix = computeStiffnessMatrix (m_elems);
	  t3 = clock();
	   cout << "Creating Stiffness and RHS consumes " << (t3-t2)/CLOCKS_PER_SEC<< endl;
	std::cout << "B =" << std::endl;
	//StiffnessMatrix.print();
	//Tinh ve phai F
	std::cout <<endl<< "F =" << std::endl;	

	F = computeRightHandSide(m_elems);	
	//F.print();
	std::cout<<"Solving C..."<<endl;
	Solver sol;
	
	R = sol.giaiHePT(StiffnessMatrix, F);
	//R.print();
	 t4 = clock();
	 cout << "Solving Linear Equation System consumes " << (t4-t3)/CLOCKS_PER_SEC<< endl;
	//R.print();
	std::cout<<ExactSolution<double>(0.356)<<" "<<app_solution(R, m_elems, 0.356)<<endl;
	cout<<"ComputeErrorAtPoint: "<<fabs(ExactSolution<double>(0.356)-app_solution(R, m_elems, 0.356))<<endl;
	t5= clock();
	 cout << "Time compute Error at point onsumes " << (t5-t4)/CLOCKS_PER_SEC<< endl;
	cout<<"ComputeL2Error: "<<computeL2Error(m_elems,R)<<endl;
	//cout<<"ComputeL22Error: "<<computeL2Error(m_elems,R,m_elems)<<endl;
	t6=clock();
	 cout << "Time compute L2 Error consumes " << (t6-t5)/CLOCKS_PER_SEC << endl;
	  cout << "Total time: " << (t6-t1)/CLOCKS_PER_SEC << "  seconds..." << endl;
	return 0;

}

FullMatrix<double> computeStiffnessMatrix(vector<FElement> ip_elems)
{	
	int N = ip_elems.size();
	FullMatrix<double> v_tmp(N - 1, N - 1);
	for (int i = 0; i < N - 2; i++)
	{
		v_tmp[i][i] = Intergration(&Bii_1, ip_elems[i], 20) + Intergration(&Bii_2, ip_elems[i+1], 20);
		v_tmp[i][i+1] = Intergration(&Bii1, ip_elems[i+1], 20);
		v_tmp[i+1][i] = v_tmp[i][i+1];
	}
	v_tmp[N - 2][N - 2] = Intergration(&Bii_1, ip_elems[N - 2], 20) + Intergration(&Bii_2, ip_elems[N-1], 20);
	return v_tmp;
}
EVector<double> computeRightHandSide(vector<FElement> ip_elems)
{
	int N = ip_elems.size();
	EVector<double> v_tmp(N-1);
	v_tmp.setZero();
	for(int i = 0; i < N - 1; i++)
	{
		v_tmp[i] = Intergration(&F_1, ip_elems[i], 20) + Intergration(&F_2, ip_elems[i+1], 20);
	}
	return v_tmp;
}
double g(double ip_x)
{
	return (B - A) * ip_x + A;
}
double norm(double (*f)(double))
{
	return Intergration(f, 20);
}
double computeL2Error(vector<FElement> ip_elems, EVector<double> R, vector<FElement> m_elems)
{
	double s = 0;
	for(int i = 0; i<N; i++)
	{
		s += Intergration(&F_error, ip_elems[i], 20, R, m_elems);
	}
	return sqrt(s);
}
double computeL2Error(vector<FElement> ip_elems,EVector<double> ip_c)
{
	assert(ip_c.size()== N-1);
	init_gauss();
	 
//    double *error = new double[GaussDegree];
	double valF=0.0;
    double x0, x1,xj, h;
    double sum = 0.0, s = 0.0;
   for(int i=0;i<N;i++)
   {
	   s = 0.0;
	  x0 =  ip_elems[i].getVertexCoord(0);
      x1 = ip_elems[i].getVertexCoord(1);
	  h = ip_elems[i].getArea();
	  for(int j=1;j<=GaussDegree;j++)
		{
			xj=x0 + xi[j][GaussDegree]*(x1-x0);
			if(i==0)
			{
				valF=ip_c[i]*(xj-x0);
			}
			if(i==N-1)
			{
				valF=(ip_c[i-1]*(x1-xj));
			}
			if((i>0)&&(i<N-1))
			{
				valF=ip_c[i-1]*(x1-xj)+ip_c[i]*(xj-x0);
			}
			valF/=h;
			valF-=ExactSolution<double>(xj);
			s+=wi[j][GaussDegree] * valF*valF;
		}
	  sum+= s*h; 
   }
   return sqrt(sum);
   return sqrt(sum*ip_elems[0].getArea());
}
double computeEError(vector<FElement> ip_elems,EVector<double> ip_c)
{
	double sum=0.0;
	double valF=0.0;
	for(int i=0;i<N-1;i++)
	{
	/*valF=ip_c[i]-*/
	}
	return sum;
}